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Abstract 

A general strategy to solve the non-perturbative renormalization problem 
in lattice QCD, using finite-size techniques and numerical simulations, is de- 
scribed. As an illustration we discuss the computation of the axial current 
normalization constant, the running coupling at zero quark masses and the 
scale evolution of the renormalized axial density. The non-perturbative cal- 
culation of 0(a) correction terms (as they appear in Symanzik's improvement 
programme) is another important field of application. 
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1. Renormalization in lattice QCD is intimately connected with the contin- 
uum limit and is thus of fundamental importance. In practical calculations, 
using numerical simulations, renormalization comes in when lattice data are 
related to quantities defined in the continuum theory. The additive and mul- 
tiplicative renormalization of the quark masses in lattice QCD with Wilson 
quarks (chosen here) is an example of this, and there are very many more 
cases where renormalization constants need to be computed. 

Perturbation theory is of limited value in this connection, because in general 
one is unable to compute more than two or three terms in the perturbation 
expansion. The numerical evaluation of the truncated series is hence liable to 
scheme ambiguities and the estimation of the truncation error becomes a matter 
of subjective judgement, particularly in situations where the coupling is large. 
Even if the series appears to converge well, some doubt will always remain that 
the result differs significantly from the exact number due to non-perturbative 
effects. 

In some cases renormalization constants can be computed through numerical 
simulations. This has always been the method of choice for the additive quark 
mass renormalization. More recently various non-perturbatively defined renor- 
malized gauge couplings have been calculated in the quenched approximation 
[1-8] and different ways to extract the renormalization constants associated 
with the isovector axial current and other composite fields have been described 
[9-15]. 

In this letter we propose to combine numerical simulations with finite-size 
techniques to solve the non-perturbative renormalization problem in lattice 
QCD. The method has previously been applied to compute the running cou- 
pling in SU(2) and SU(3) gauge theories [1-5]. Our aim here is to advertise 
its use in QCD by explaining the basic strategies in the context of this theory 
and by going through a list of applications. 

2. When computing renormalization constants through numerical simula- 
tions, one is confronted with a number of technical difficulties. One of these 
arises from the fact that renormalization is scale dependent in general. It is 
then often necessary to trace the evolution of the renormalized parameters and 
renormalized composite fields from low to high energies, where contact can 
be made with perturbative renormalization schemes. For lattice gauge theory 
this presents a problem, since one cannot afford to simulate lattices covering 
physical scales orders of magnitude apart. 

A second difficulty derives from the well-known limitation that simulations 
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of QCD with very light quarks are usually prohibitively costly. As a result 
one is forced to adopt a mass dependent renormalization scheme or to rely on 
extrapolations to the chiral limit. Both options are not particularly attractive 
and it would be much preferable, if one were able to perform the parameter 
and field renormalizations directly at vanishing quark masses. 

The infiuence of lattice effects is a further problem that should be addressed 
when discussing non-perturbative renormalization. A good example to illus- 
trate this is the renormalization constant Z\ of the (local) isospin current 
[9-14]. Z\ may be fixed by requiring the isospin charge of some low-energy 
states to assume integer values. Depending on exactly which states are chosen 
to compute Z\ ^ results differing by terms of order a (where a denotes the lat- 
tice spacing) are obtained. Such ambiguities may be considered a systematic 
error on Z\ ^ but this point of view is not entirely satisfactory, since it is dif- 
ficult to estimate the error in an objective manner. A better way to deal with 
the problem is to fix Z\ through a definite normalization prescription and to 
study the approach to the continuum limit of the physical quantities involving 
the renormalized vector current that one is interested in. 

Moving closer to the continuum limit means larger lattices and hence rapidly 
increasing cost. Less expensive ways to reduce cutoff effects include the use of 
an 0(a) improved action [16-20] along with improved composite fields [21]. A 
problem which one has here is that improvement needs to be verified. Moreover 
one may not be satisfied with a perturbative computation of the coefficients 
multiplying the 0(a) correction terms added to the action and the composite 
fields. For optimal improvement a non-perturbative determination of these 
coefficients may be necessary. 

The finite-size technique discussed below is able to cope with all three prob- 
lems mentioned in this section. This is rather non-trivial and we shall not 
attempt to present the method in full detail here. Instead we shall proceed 
from simple to more complex applications and defer all technical discussions 
to later publications. 

3. We begin by considering QCD in a euclidean space-time volume of size 
T X , where the time-like extent T and the spatial size L are taken to be in 
a fixed proportion {T = 2L, for example). With suitable boundary conditions 
the renormalization of the theory in such a world proceeds as in ordinary space- 
time. In particular, the counterterms that must be added to the bare action 
may be chosen to be independent of L and we shall assume that this is what 
has been done. From the point of view of the lattice theory, this simply means 
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that the continuum limit at fixed L (given in physical units) is obtained with 
the usual scaling of the bare coupling and quark masses. 

A finite volume is usually regarded as a source of systematic error in nu- 
merical simulations of lattice QCD. Here, however, it is taken as a probe of 
the theory. To study renormalization and Symanzik improvement the relevant 
sizes L are then smaller than 1 fm (hence the name "femto-universe" [22]). We 
shall in fact take L to very small values when tracing the scale evolution of the 
renormalized parameters and fields to high energies. 

Evidently the femto-universe is experimentally inaccessible and one may, 
therefore, be driven to conclude that it can only be of academic interest. This 
argumentation overlooks the fact that there are theoretical ways to relate finite 
volume with infinite volume physics, either through renormalized perturbation 
theory (at high energies) or via the bare parameters and fields on the lattice. 
It is then possible to arrange the computations so that all reference to a finite 
volume disappears from the final results. In other words, the femto-universe is 
only used as an intermediate device. 

So far we did not specify the boundary conditions on the quark and gauge 
fields. Different boundary conditions amount to probing the theory in different 
ways, i.e. there are no fundamental reasons for choosing any particular pre- 
scription. To avoid inessential technical complications it is however wise to 
require that zero modes are excluded at tree-level of perturbation theory. A 
possible choice then are periodic boundary conditions in all spatial directions 
and Dirichlet boundary conditions at time xq = Q and xq = T . Explicitly, at 
xq = Q the gauge potential A^{x) is required to satisfy f 

Ak{x) = Ck{-^), A; = 1,2, 3, (3.1) 

where Ck is a given classical field. Since the Dirac equation is of first order, 
only half of the components of the quark and anti-quark fields il}{x) and il}{x) 
can be prescribed, viz. 

i(l + 7o)V'(a;) = />(x), V^(a;)(l-7o)i = />(x). (3.2) 

Similar boundary conditions are imposed sX xq = T [23]. 

The euclidean functional integral, which may now be set up in the usual way, 
depends on the boundary values of the fields and may be interpreted as the 

\ For ease of presentation we use a continuum notation in this paper, but the formulae 
should be taken as shorthand for the corresponding lattice expressions 
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quantum mechanical amplitude for going from the given field values at time 
a;o = to the values specified at time xq = T. It is, therefore, also referred to 
as the "Schrodinger functional". The renormalizability of the theory in these 
conditions has been studied in refs. [2,23,24] and it was found that no new coun- 
terterms are needed except for one term at the boundary which amounts to 
rescaling the boundary values p and by a logarithmically divergent renor- 
malization constant. The renormalization of the coupling, the quark masses 
and the composite fields at times < xq < T are not affected by the presence 
of the boundaries. 

One of our motivations for choosing boundary conditions as specified above 
was to ensure that no zero modes occur to lowest order of perturbation the- 
ory. The eigenvalues of the free Dirac operator are in fact separated from zero 
by a gap which is proportional to the infrared cutoff 1/L at vanishing quark 
masses [23]. We have verified through numerical simulations, using a recently 
developed eigenvalue program [25], that the gap persists in the fully interact- 
ing theory for box sizes L < 1 fm. Moreover the fiuctuations of the lowest 
eigenvalues remain small compared to the gap. The important consequence of 
this observation is that numerical simulations of the Schrodinger functional are 
feasible at zero quark masses. In the interesting range of L such simulations 
are not significantly more expensive than simulations at small positive quark 
masses. 

As an aside we remark that the Dirac operator in the continuum theory, in 
an arbitrary smooth background gauge field and with Schrodinger functional 
boundary conditions, can be shown to have no zero modes at vanishing quark 
masses. In particular, the Atiyah-Singer index theorem does not apply when 
such boundary conditions are imposed. 

4. In the standard formulation of lattice QCD with two fiavours of mass- 
degenerate Wilson quarks there are two bare parameters, the gauge coupling 
^0 and the quark mass mo. The renormalized quark mass nifi is related to the 
bare mass through 



where and nic are the multiplicative and additive mass renormalization 
constants. We now describe how to compute nic using the femto-universe. 
The idea is to consider the PCAC relation 
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between the isovector axial current 



A;{x) = i;{x)j,j5^T-i;{x) (4.3) 

and the associated (unrenormalized) density 

P^x) = ^{x)j,^T^i,{x) (4.4) 

(r" denotes a Pauli matrix acting on the flavour indices of the quark flelds). 
The mass m appearing in eq.(4.2) is proportional to the renormalized quark 
mass. In particular, the critical bare mass nic may be calculated by flnding 
the value of nio where m vanishes. 

For given bare parameters the mass m can be extracted from ratios of cor- 
relation functions, 

m = i {d,Al{x)0^) I {P-{x)0^) + 0(a), (4.5) 

where (9" is a suitable polynomial in the quark and gluon flelds supported in 
a region of space-time not containing x. Up to cutoff effects of order a, the 
computed value of m should of course be independent of the source (9" and 
also of the volume and the boundary values of the flelds |. 

An interesting set of sources (9" may be constructed by differentiating the 
Schrodinger functional with respect to the boundary values of the quark held. 
This amounts to inserting quark flelds at the boundary and so we introduce 
symbolic boundary flelds C(x) ^-nd C(x) through 

In eq.(4.5) we now choose 

0^= rdVd3zC(y)75KC(z) (4.7) 
Jo 

and set the boundary values p and p to zero after differentiation. The boundary 
values of the gauge held are taken to be zero or constant abelian as in ref.[5]. 



I Recall that eq.(4.2) is a special case of the euclidean field equations. These are derived 
locally from the action and so do not depend on the boimdary conditions (as long as one 
keeps away from the boundaries) 
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Fig. 1. Values of m as computed from eq.(4.5), using tlie Wilson quark action. 
The open and full symbols correspond to zero and non-zero boundary values of 
the gauge field. 

Note that the boundary fields in eq.(4.7) are projected to their zero momentum 
components. We do not need to choose a gauge to do this, since the gauge 
symmetry at the boundary is fixed by the boundary value of the gauge field. 

The computation of the correlation functions in eq.(4.5) through numerical 
simulation is straightforward. A typical result from a 16 X 8^ lattice with 
quenched Wilson quarks is shown in fig. 1. The bare mass nio is close to nic in 
this example and the bare coupling is such that d/g^ = 6.4, which corresponds 
to a physical lattice size L of about 0.4 fm [8,3]. Naively one would expect 
to see a plateau in this plot, since m should be independent of the time xq at 
which the axial current is inserted. This is far from being the case. Moreover if 
we change the boundary values of the gauge field from zero to some weak field 
(half the strength of the field used in the computation of the running coupling 
in ref.[3]), the value of m in the middle of the lattice is shifted by as much as 
32 MeV. 

Further studies reveal that these unexpected findings result from violations 
of chiral symmetry by lattice effects. The strongest effects (those of order a) 
can be canceled by using on-shell improved fields and an improved action. The 
latter is obtained by adding a Pauli term, 

^Cs„a(7^^F^^(a;), (4.8) 
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Fig. 2. Same as fig. 1, but using tlie 0(a) improved action and tlie improved 
axial current A1^/Za and density / Zp witli Cgw = 1-60 and ca = —0.027. 

to the lattice Dirac operator, where F^ij{x) denotes the gluon field strength 
[19,20]. For the renormalized on-shell 0(a) improved axial current and density 
we may take 

Al = Z^{{l + h^am^)Al + c^ad^P°'] , m^ = mQ-m,, (4.9) 
p« = Zp(l + 5pamq)P". (4.10) 



With these modifications PCAC is expected to be valid up to errors of order 
a^, provided the improvement coefficients Cs^v, 5a j and ca are chosen appro- 
priately. To lowest order of perturbation theory Cs^v = 5a = = 1 and ca = 
[19,21]. 

The effect of improvement is quite dramatic (see fig. 2). There is now a 
wide plateau and the remaining dependence on the boundary values of the 
gauge field is barely significant. The 0(a) corrections associated with 5a and 
5p are unimportant in the present context since they just amount to a rescaling 
of m by a constant factor. They are anyway numerically insignificant at small 
quark masses and we have thus decided to neglect them in this calculation. 

The quoted values of Cs^v and ca have been obtained by computing m from 
three different correlation functions and adjusting the improvement coefficients 
until the same result is obtained in all three cases. In other words, the coef- 
ficients are determined by imposing a non-perturbative "improvement condi- 
tion", which amounts to requiring PCAC to hold exactly in three correlation 
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Fig. 3. Computation of tlie critical bare mass nic in tlie 0(a) improved 
tlieory. Tlie lattice parameters are the same as in fig. 2. The three data points 
are interpolated linearly and nic is found at m = 0. 

functions, with the same value of m. It should again be emphasized at this 
point that the improved action and fields are independent of the physical situa- 
tion considered. Our results for the improvement coefficients are hence directly 
relevant for the computation of low-energy properties of QCD in large volumes. 

After improvement has been taken into account the calculation of the critical 
bare mass nic proceeds along the lines described earlier in this section. From 
the plateau in fig. 2 the mass m is obtained with small errors. Repeating this 
calculation for two more values of the bare mass then yields the data points 
shown in fig. 3 and nic is finally obtained by linear interpolation. If there 
should be any doubt about the reliability of the interpolation, one can always 
verify that m = within statistical errors at the computed value of nic- The 
slope of the line in fig. 3 is about 1.09 and ni is hence close to nio — nic in the 
small quark mass region. 

5. The axial current normalization constant Z\ appearing in eq.(4.9) is 
determined by the chiral Ward identities. This has previously been discussed 
in a series of publications [9-13] and practical methods to extract Z\ from the 
identities have been described, using numerical simulations. 

The chiral Ward identities are derived locally from the QCD action and so 
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are valid in finite volumes, too. The computation of Zp^ can thus be carried 
out in the femto-universe in essentially the same way as on physically large 
lattices. If we set the quark masses to zero (which simplifies the calculation 
considerably), the relevant Ward identity may be written in the form 

/ Aa,{x) e»^= {Al{x)A\{y)Q'^) = 2^ (V=(y)Q=) + 0{a'). (5.1) 

■JdR 

The integral in this formula is taken over the boundary of some space-time 
region R containing the point y. is an arbitrary source located outside R 
and V^(a;) denotes the renormalized on-shell improved isospin current, 

V;{x) = Zy{{l + 5vamq)y; + cyad^T;^} , (5.2) 

V;{x)=lp{x)j^^T-^{x), T;,{x)=lp{x)a^,^T-^{x). (5.3) 

Note that the 0(a) tensor correction to the vector current does not contribute to 
the isospin charge, but it could be important when computing fp, for example. 

A convenient choice for R is the region between two equal time hyper-planes. 
For the source one may take 

where O"" is given by eq.(4.7) and O'"" is constructed similarly at xq = T . If 
we set z/ = in eq.(5.1), the correlation function on the right hand side reduces 
to a matrix element of the isospin charge between states with isospin 1 and so 
is determined by the normalization condition for the vector current. On the 
other side of the equation the normalization constant Zp^ is the only unknown 
coefficient, i.e. it can be computed by requiring eq.(5.1) to hold exactly for the 
chosen source and lattice parameters. 

An attractive feature of this computation is that the fields in the Ward iden- 
tity (5.1) are localized in different space-time regions. Cutoff effects from 
contact terms are hence avoided and the error can be argued to be of order , 
if only the on-shell improved currents and the improved action are employed. 

6. So far we have been exclusively concerned with scale independent renor- 
malization constants. In this last section we address the more difficult problem 
of scale dependent renormalization. It is here that the proposed finite-size tech- 
nique develops its full power. 
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Conventionally the renormalization of QCD is discussed in perturbation the- 
ory. In this framework the currently most popular renormalization scheme is 
the MS scheme of dimensional regularization. The associated renormalization 
scale 1^ is introduced in a rather implicit manner, but it is clear that it should 
be large compared to (say) the nucleon mass so that perturbation theory may 
be expected to apply. 

At low energies, on the other hand, it is natural (and common practice 
in lattice QCD) to fix the bare parameters in the action by requiring that the 
nucleon mass and the masses of some pseudo-scalar mesons (one for each quark 
fiavour) assume prescribed values. It is also possible to define renormalized 
composite fields by specifying some of their hadronic matrix elements. Such 
renormalization schemes are referred to as hadronic schemes. 

An important part of the non-perturbative renormalization problem in lattice 
QCD is to determine the relation between the hadronic and the perturbative 
schemes. In particular, one would like to be able to compute the running 
coupling and the running quark masses at high energies given the nucleon 
and meson masses. One is also often confronted with the task of computing 
hadronic matrix elements of some local composite fields whose normalization 
is given at high energies through the MS scheme. 

As already mentioned in sect. 2, the principal difficulty in such computations 
is to relate physical quantities defined at energy scales differing by factors of 
10 to 100 (and may be more). In the following we describe in outline how this 
problem can be solved using finite-size techniques. The method has first been 
discussed in ref.[l] and was later applied in pure gauge theories [2-5]. 

The basic idea is to introduce an intermediate finite-volume renormaliza- 
tion scheme, where all renormalized quantities are defined at scale 1/L. The 
Schrodinger functional provides a technically attractive framework to set up 
such a renormalization scheme, which is then called the SF scheme. Using 
perturbation theory the SF scheme can be matched with the MS scheme when 
L is sufficiently small and /^L of order 1. At the lower end of the energy scale, 
the SF scheme can be related to the hadronic schemes through numerical sim- 
ulations. The crucial point now is that the scale evolution in the SF scheme 
is computable through a non-perturbative recursive procedure. Via the SF 
scheme one is hence able to connect the MS scheme with the hadronic schemes. 

We now need to be a bit more explicit on the definition of the SF scheme 
and shall then explain how the scale evolution of the renormalized parameters 
and fields can be computed. We choose to set up the SF scheme in such a way 
that all normalization conditions are given at zero quark mass [26]. 
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The renormalized coupling aspig), q = 1/-^, is obtained from the Schro- 
dinger functional by calculating the response of the system to changes of the 
boundary values of the gauge field (full details are given in refs.[3,24]). We set 
rriQ = rric in this definition, as indicated above, and scale the chosen boundary 
values proportionally to L. The only variables on which asF(?) depends are 
then the bare coupling and the lattice size L/a. It has been shown that 
«SF (?) is a renormalized quantity which satisfies the expected massless renor- 
malization group equation. Moreover its relation to Oiy^{q) has been worked 
out to one-loop order [24]. 

For the running quark mass we may take 

msF{q) = mZx/Zp, g = (6.1) 

where m is determined from PCAC using the (unrenormalized) improved cur- 
rent A^^/Zx and density V"" jZ^. The computation of the axial current nor- 
malization constant Zp^ has been discussed in the previous section. To fix the 
renormalization constant Z^ of the axial density (4.10), we must impose a 
normalization condition. A simple possibility is to require that 

(P"(a;)C")2 = -(9/L6)(C)'"C"), (6.2) 

at zero quark mass, zero boundary values, = r/2 and T = 2L. The 
proportionality factor in this equation has been chosen so that = 1 to lowest 
order of perturbation theory. Other composite fields can be renormalized in a 
similar manner. 

The evolution of asF(?) from some large initial value of q towards lower 
scales may now be calculated simply by increasing L at fixed bare parameters. 
Of course we can only increase L by a factor of 2 or so, as otherwise one ends 
up with lattices that cannot be simulated with the currently available comput- 
ers. In a second step the lattice spacing a is increased at fixed renormalized 
parameters. This amounts to lowering Lj a by say a factor of 2 and adjusting 
the bare parameters so that the quark mass (equal to zero) and the renormal- 
ized coupling remain unchanged. Now we can evolve the coupling by another 
factor of 2, and so on, until the low-energy regime is reached. 

With little additional effort one can also compute the scale evolution of the 
renormalization constant Z^ (which in turn determines the evolution of the 
running quark mass). In the recursion we simply have to calculate the ratios 
Zp{go, L'/a)/Zp{go, L/a) when the box size is increased from L to L' at fixed 
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bare parameters. The product of all these ratios is equal to the total change 
in the normalization of the renormalized axial density during the evolution. It 
should be emphasized that the lattice cutoff 1 j a is always much greater than the 
renormalization scale, at all stages of the calculation. Cutoff effects are hence 
expected to be small. Moreover, following refs.[l-5], any remaining effects 
can be extrapolated away by simulating sequences of lattices with decreasing 
lattice spacings. 

Before the running coupling and quark mass can be computed along the lines 
explained above, some preparatory work is still required (simulation algorithms 
for dynamical quarks with 0(a) improved action need to be developed, for ex- 
ample). As an intermediate step one may be interested to perform simulations 
with zero and negative numbers of dynamical quark flavours [27,28], since these 
are relatively cheap and thus provide a useful laboratory to study some of the 
technical issues involved. Such calculations may also give interesting additional 
insights into the dependence of the evolution of the renormalized parameters 
and fields on the number of light quarks. 

We are indebted to Roberto Frezzotti and Chris Sachrajda for clarifying dis- 
cussions on the improved axial current. The numerical simulations have been 
performed on the powerful APE/Quadrics computers at DESY-IfH (Zeuthen). 
We thank the staff of the computer center at Zeuthen for their support. 
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